Curso II – SISSA

Fusión de Productos de Precipitación con datos in-situ

Objetivo y contenido

Objetivo y contenido

Objetivo: Fusionar productos de precipitación con datos de estaciones

  • ¿Por qué fusionar productos de precipitación y datos de estaciones?
  • Consideraciones a tener en cuenta
  • Metodología de fusión
  • Ejemplo de fusión de productos de precipitación con datos de estaciones

Fusión de Productos de Precipitación

Fusión de Productos de Precipitación

La fusión de productos de precipitación con datos in-situ tiene como objetivo combinar y aprovechar la información proporcionada por ambos tipos de datos para mejorar la calidad y la resolución espacial de las estimaciones de precipitación.

Existen varias razones por las cuales se fusionan los productos de precipitación con los datos in-situ:

  • Complementariedad de datos
  • Mejora de la resolución espacial

Fusión de Productos de Precipitación

En el ejemplo de fusión de productos de precipitación con datos in-situ, se aplicará una metodología específica para combinar los datos y se mostrará cómo se pueden obtener estimaciones de precipitación mejoradas utilizando ambos tipos de datos.

Consideraciones

Consideraciones

La fusión de productos de precipitación con datos in-situ requiere considerar varios puntos:

  • homogeneidad de las fuentes de datos,

  • calidad y confiabilidad de los datos in-situ

Se pueden utilizar diferentes enfoques de fusión, como métodos estadísticos, interpolación espacial o modelos de fusión.

Metodología

Baez-Villanueva et al., (2020)

Ejemplo de fusión de productos de precipitación con datos de estaciones

Ejemplo

En este ejemplo, fusionaremos dos productos de precipitación (CHIRPSv2 y PERSIANN-CDR) con estaciones in situ en Valparaíso, Chile.

Paso 1: Cargar los paquetes necesarios para trabajar con las estaciones:

library(zoo)
library(sf)
library(rgdal)
library(raster)

El paquete RFmerge se puede instalar desde Github:

remotes::install_github('hzambran/RFmerge')
library(RFmerge)

Paso 2: Cargar los datos de las estaciones y el shapefile del área de estudio. Para esto, almacenaremos las rutas de los archivos:

observations <- 'C:/Users/Data/Timeseries.zoo'
metadata     <- 'C:/Users/Data/Metadata.csv'
shape        <- 'C:/Users/Data/Valparaiso.shp'

En este paso, estamos especificando las rutas de los archivos de datos de las estaciones y del shapefile del área de estudio. Asegúrate de proporcionar las rutas correctas a tus archivos correspondientes.

Podemos utilizar la función read.zoo del paquete zoo para importar las series temporales, la función read.csv para importar los metadatos y la función read_sf del paquete sf para importar el shapefile.

observations <- read.zoo(observations, header = TRUE)
metadata     <- read.csv(metadata)
shape        <- read_sf(shape)

Paso 3: Ahora podemos cargar los datos diarios de los productos de precipitación y el modelo digital de elevación. Para esto, almacenaremos las rutas de los archivos correspondientes.

chirps        <- 'C:/Users/Data/P_products/CHIRPS5km.tif'
persiann_cdr  <- 'C:/Users/Data/P_products/PERSIANNcdr5km.tif'
dem           <- 'C:/Users/Data/DEM/ValparaisoDEM5km.tif'

Podemos utilizar la función brick del paquete raster para importar los productos de precipitación (archivos multibanda) y la función raster para importar el DEM.

chirps        <- brick(chirps)
persiann_cdr  <- brick(persiann_cdr)
dem           <- brick(dem)
print(chirps)
## class      : RasterBrick 
## dimensions : 40, 38, 1520, 243  (nrow, ncol, ncell, nlayers)
## resolution : 0.05, 0.05  (x, y)
## extent     : -71.85, -69.95, -34, -32  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : CHIRPS5km.tif 
## names      : CHIRPS5km_1, CHIRPS5km_2, CHIRPS5km_3, CHIRPS5km_4, CHIRPS5km_5, CHIRPS5km_6, CHIRPS5km_7, CHIRPS5km_8, CHIRPS5km_9, CHIRPS5km_10, CHIRPS5km_11, CHIRPS5km_12, CHIRPS5km_13, CHIRPS5km_14, CHIRPS5km_15, ...
print(persiann_cdr)
## class      : RasterBrick 
## dimensions : 40, 38, 1520, 243  (nrow, ncol, ncell, nlayers)
## resolution : 0.05, 0.05  (x, y)
## extent     : -71.85, -69.95, -34, -32  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : PERSIANNcdr5km.tif 
## names      : PERSIANNcdr5km_1, PERSIANNcdr5km_2, PERSIANNcdr5km_3, PERSIANNcdr5km_4, PERSIANNcdr5km_5, PERSIANNcdr5km_6, PERSIANNcdr5km_7, PERSIANNcdr5km_8, PERSIANNcdr5km_9, PERSIANNcdr5km_10, PERSIANNcdr5km_11, PERSIANNcdr5km_12, PERSIANNcdr5km_13, PERSIANNcdr5km_14, PERSIANNcdr5km_15, ...
print(dem)
## class      : RasterBrick 
## dimensions : 40, 38, 1520, 1  (nrow, ncol, ncell, nlayers)
## resolution : 0.05, 0.05  (x, y)
## extent     : -71.85, -69.95, -34, -32  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : ValparaisoDEM5km.tif 
## names      : ValparaisoDEM5km 
## min values :          11.4891 
## max values :         5123.841

Los dos archivos geotiff multibanda proporcionados no almacenan la fecha de la estimación de precipitación. Por lo tanto, antes de realizar cualquier análisis exploratorio, es necesario asignar nombres significativos a cada banda. Esto nos permitirá tener una mejor idea de las fechas asociadas a cada capa y facilitará el análisis posterior de los datos.

ldates <- hydroTSM::dip("1983-01-01", "1983-08-31")
## Please note that 'maptools' will be retired during October 2023,
## plan transition at your earliest convenience (see
## https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
## for guidance);some functionality will be moved to 'sp'.
##  Checking rgeos availability: FALSE
names(chirps)       <- ldates
names(persiann_cdr) <- ldates

A continuación, queremos visualizar las primeras y últimas seis filas de los metadatos espaciales:

head(metadata)
##       Code      lon      lat            NOM_REG  NOM_PROV COD_COM     NOM_COM
## 1 P5101005 -70.8000 -32.0836 Regin de Valparaso Los Andes    5304 San Esteban
## 2 P5111002 -71.0300 -32.1567 Regin de Valparaso Los Andes    5304 San Esteban
## 3 P5101006 -70.7833 -32.1836 Regin de Valparaso Los Andes    5304 San Esteban
## 4 P5100006 -70.7839 -32.2261 Regin de Valparaso Los Andes    5304 San Esteban
## 5 P5100005 -70.7100 -32.2286 Regin de Valparaso Los Andes    5304 San Esteban
## 6 P5111001 -71.1386 -32.2528 Regin de Valparaso Los Andes    5304 San Esteban
tail(metadata)
##        Code      lon      lat            NOM_REG  NOM_PROV COD_COM     NOM_COM
## 29 P5427006 -71.2144 -33.0986 Regin de Valparaso Los Andes    5304 San Esteban
## 30 P5510002 -71.5553 -33.1450 Regin de Valparaso Los Andes    5304 San Esteban
## 31 P5741002 -71.1467 -33.1686 Regin de Valparaso Los Andes    5304 San Esteban
## 32 P5530002 -71.6250 -33.5747 Regin de Valparaso Los Andes    5304 San Esteban
## 33 P5748003 -71.5106 -33.6344 Regin de Valparaso Los Andes    5304 San Esteban
## 34  P330030 -71.6142 -33.6550 Regin de Valparaso Los Andes    5304 San Esteban

Ahora podemos graficar la serie temporal diaria de precipitación para la primera estación (código: P5101005).

main <- paste("Precipitación diaria para la estación", metadata$Code[1])
ylab <- "Precipitación [mm]"
x_ts <- observations[,1]
#plot(x.ts, main=main, ylab= ylab, col="blue")
hydroTSM::hydroplot(x_ts, pfreq="o", main=main, ylab= ylab)
## [Note: pfreq='o' => ptype has been changed to 'ts']
grid()

Podemos graficar las estimaciones de precipitación acumulada para los primeros ocho meses de 1983 de CHIRPS y PERSIANN-CDR, y superponer los límites del área de estudio:

chirps_total   <- sum(chirps, na.rm= FALSE)
persiann_total <- sum(persiann_cdr, na.rm= FALSE)
plot(chirps_total, main = "CHIRPSv2 [Enero - Agosto] ", xlab = "Longitud", ylab = "Latitud")
lines(shape)

plot(persiann_total, main = "PERSIANN-CDR [Enero - Agosto] ", xlab = "Longitud", ylab = "Latitud")
lines(shape)

Paso 4: Convertir el objeto metadata en un objeto espacial. Para esto, lo transformatemos en un SpatialPointsDataFrame utilizando los campos de latitud y longitud, almacenados en las columnas lat y lon:

stations <- st_as_sf(metadata, coords = c('lon', 'lat'), crs = 4326)
print(stations)
## Simple feature collection with 34 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -71.625 ymin: -33.655 xmax: -70.4719 ymax: -32.0836
## Geodetic CRS:  WGS 84
## First 10 features:
##        Code            NOM_REG  NOM_PROV COD_COM     NOM_COM
## 1  P5101005 Regin de Valparaso Los Andes    5304 San Esteban
## 2  P5111002 Regin de Valparaso Los Andes    5304 San Esteban
## 3  P5101006 Regin de Valparaso Los Andes    5304 San Esteban
## 4  P5100006 Regin de Valparaso Los Andes    5304 San Esteban
## 5  P5100005 Regin de Valparaso Los Andes    5304 San Esteban
## 6  P5111001 Regin de Valparaso Los Andes    5304 San Esteban
## 7  P5110003 Regin de Valparaso Los Andes    5304 San Esteban
## 8  P5111004 Regin de Valparaso Los Andes    5304 San Esteban
## 9  P5120004 Regin de Valparaso Los Andes    5304 San Esteban
## 10 P5200006 Regin de Valparaso Los Andes    5304 San Esteban
##                     geometry
## 1     POINT (-70.8 -32.0836)
## 2    POINT (-71.03 -32.1567)
## 3  POINT (-70.7833 -32.1836)
## 4  POINT (-70.7839 -32.2261)
## 5    POINT (-70.71 -32.2286)
## 6  POINT (-71.1386 -32.2528)
## 7  POINT (-70.9997 -32.2756)
## 8  POINT (-71.0792 -32.3058)
## 9  POINT (-71.2425 -32.3075)
## 10 POINT (-70.7528 -32.3408)

Ahora podemos graficar el DEM con las ubicaciones de las estaciones pluviométricas que se utilizarán en el cálculo de RF-MEP, y superponiendo los límites del área de estudio y las estaciones:

plot(dem, main="SRTM-v4", xlab="Longitud", ylab="Latitud", col=terrain.colors(255))
lines(shape)
plot(st_geometry(stations), pch = 16, col="black", add = TRUE)

Paso 5: Queremos producir un producto fusionado de precipitación desde el 1 de enero hasta el 31 de agosto de 1983 (243 días). Por lo tanto revisaremos la longitud de los productos.

nlayers(chirps)
## [1] 243
nlayers(chirps) == nlayers(persiann_cdr)
## [1] TRUE
nlayers(chirps) == nrow(observations)
## [1] TRUE

Además, revisaremos la extensión espacial y la resolución espacial de los productos.

extent(chirps) == extent(persiann_cdr)
## [1] TRUE
res(chirps)    == res(persiann_cdr)
## [1] TRUE TRUE

Paso 6: Ahora vamos a fusionar los productos de precipitación, el DEM y los datos de las estaciones. Para esto, prepararemos una lista de covariables.

covariates <- list(chirps=chirps, persiann_cdr=persiann_cdr, dem=dem)

Si deseas que los archivos resultantes fusionados se escriban en disco, debes definir el directorio de salida (drty.out) antes de ejecutar RFmerge. Luego, puedes ejecutar la función RFmerge de la siguiente manera:

drty_out <- file.path("C:/Data/Merging/Output/")
rfmep <- RFmerge(x=observations, metadata=metadata, cov=covariates,
                 id="Code", lat="lat", lon="lon", mask=shape,
                 training=0.8, write2disk=TRUE, drty.out=drty_out)

Si write2disk=TRUE, los siguientes archivos de salida se almacenarán en el directorio drty.out definido por el usuario:

    1. las estaciones de medición de precipitación utilizadas como conjuntos de datos de entrenamiento y evaluación; y
    1. el producto fusionado final (archivos GeoTiff individuales).

Los objetos resultantes mencionados anteriormente se almacenarán dentro de drty.out de la siguiente manera:

  • Ground_based_data/Training/: datos de series temporales y metadatos utilizados para entrenar RF-MEP como un archivo zoo (Training_ts.txt) y un archivo de texto (Training_metadata.txt), respectivamente.

Los objetos resultantes mencionados anteriormente se almacenarán dentro de drty.out de la siguiente manera:

  • Ground_based_data/Evaluation/: Este directorio almacenará las series temporales y metadatos que no se utilizaron para entrenar el modelo Random Forest, pero que están disponibles para evaluar el rendimiento de los resultados en un conjunto de datos de evaluación independiente. Los archivos Evaluation_ts.txt y Evaluation_metadata.txt se almacenan como archivos zoo y CSV, respectivamente.

Los objetos resultantes mencionados anteriormente se almacenarán dentro de drty.out de la siguiente manera:

  • RF-MEP/: Este directorio almacenará los archivos GeoTiff individuales producidos por el algoritmo RF-MEP, utilizando la misma resolución espacial que las covariables seleccionadas.

rfmep_total <- sum(rfmep, na.rm = TRUE)
totals <- stack(chirps_total, persiann_total, rfmep_total)
names(totals) <- c('CHIRPSv2', 'PERSIANN_CDR', 'RF-MEP')
plot(totals)

Paso 7: utilizaremos el conjunto de datos de evaluación de observaciones de estaciones de medición de precipitación para evaluar el rendimiento de RF-MEP con observaciones que no se utilizaron en el procedimiento de fusión.

eval          <- paste0(drty_out, "/Ground_based_data/Evaluation/Evaluation_ts.txt")
metadata_eval <- paste0(drty_out, "/Ground_based_data/Evaluation/Evaluation_metadata.txt")
eval          <- read.zoo(eval, header = TRUE)
metadata_eval <- read.csv(metadata_eval)
stations_eval <- st_as_sf(metadata_eval, coords = c('lon', 'lat'), crs = 4326)

Visualicemos una capa de RF-MEP junto con las estaciones de medición de precipitación y el área de estudio.

plot(rfmep[[11]], main = "Precipitación RF-MEP para 1983-01-11", xlab="Longitude", ylab="Latitude")
lines(shape)
plot(st_geometry(stations_eval), pch = 16, col="black", add = TRUE)

Primero, compararemos RF-MEP (y los dos productos de precipitación utilizados como covariables) con los datos de las estaciones de lluvia del conjunto de evaluación. Para extraer los valores de precipitación de RF-MEP en las ubicaciones de las estaciones, utilizaremos la función extract del paquete raster:

rfmep_ts    <- t(raster::extract(rfmep, stations_eval))
chirps_ts   <- t(raster::extract(chirps, stations_eval))
persiann_ts <- t(raster::extract(persiann_cdr, stations_eval))

Ahora realizaremos una evaluación similar a la que hicimos en el Módulo 6 para evaluar el rendimiento de los productos de precipitación:

products  <- list(chirps_ts, persiann_ts, rfmep_ts)
nproducts <- length(products)
nstations <- ncol(eval)
tmp       <- rep(NA, nstations)
kges      <- data.frame(ID=stations_eval[["Code"]], CHIRPS=tmp, PERSIANN_CDR=tmp, RF_MEP=tmp)

for (i in 1:nproducts) {
  ldates        <- time(eval)
  lsim          <- zoo(products[[i]], ldates)
  kges[, (i+1)] <- hydroGOF::KGE(sim= lsim, obs= eval, method="2012")
}

Ahora podemos visualizar los resultados:

print(kges)
##         ID    CHIRPS PERSIANN_CDR    RF_MEP
## 1 P5120004 0.2456329    0.2568213 0.7087269
## 2 P5120003 0.4103445    0.3358265 0.8891823
## 3 P5414004 0.4492111    0.2953030 0.9106920
## 4 P5410007 0.4403373    0.1509242 0.9001426
## 5 P5423013 0.2988270    0.3147489 0.8865164
## 6 P5426004 0.3098025    0.3046178 0.5345708
## 7 P5427007 0.2280326    0.2647549 0.8525314

Finalmente, se generará un boxplot que compare el desempeño, en términos de KGE’, del producto fusionado en comparación con los productos originales utilizados como covariables.

sres_cols <- c("powderblue", "palegoldenrod", "mediumseagreen")
boxplot(kges[,2:4], main = "Evaluación de KGE para enero - agosto de 1983",
xlab = "Productos de precipitación", ylab = "KGE'", ylim = c(0, 1), col=sres_cols)
grid()
legend("topleft", legend=c("CHIRPS", "PERSIANN-CDR", "RF-MEP"), col=sres_cols,
pch=15, cex=1.5, bty="n")
grid()

¡Gracias por su atención!